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Jin et al reported that axisymmetric simulations of NS-like objects with polytropic EOS undergo 
critical gravitational collapse. As the critical collapse observed via fine-tuning of the adiabatic index 
F, they conjecture that critical phenomena may occur in realistic astrophysical scenarios. To clarify 
the implications this numerical observation has on realistic astrophysical scenarios, here, we perform 
dynamical analysis on the structure of the critical collapse observed in the former work. We report 
the time scales and oscillation frequencies exhibited by the critical solution and compare these results 
with values obtained from analytic perturbative mode analysis of equilibrium TOV configurations. 
We also establish the universality of the critical solution with respect to a 1-parameter family of 
initial data as well as the phase space manifold of the critical collapse. 
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,— — i Seel. Introduction. Critical phenomena in gravitational collapse was discovered by Choptuik almost a decade 

ago pQ. Since then, there has been much development in establishing the theory in its mathematical aspects. Past 
studies have considered the assumption of both spherically symmetric and axisymmetric systems, and several matter 
models eg neutron stars, perfect fluids, real scalar fields, 2D Sigma models, SU(2) Yang-Mills and Skyrme models, 
were considered in these studies. However, thus far, none have considered axisymmetric systems of neutron stars. 

i— i The neutron star models that have been considered [5] employ spherical symmetry and Type II critical phenomena 
jZ-i were found in these studies. 

Jin et al [2] carried out axisymmetric simulations of neutron stars and found Type I critical phenomena ie. black 
tTp. holes with a mass gap form in the black hole phase. The axisymmetry in these studies enable high-resolution finite 

i__i differencing that contributes to the observation of the fine structure of critical phenomena when n-paramctcrs of the 
system are fine-tuned. In this study, Jin et al tuned the central densities, the boost velocities, the initial separations 
and the polytropic index of the EOS of the NSs one at a time. Varying solely the central densities of the NSs 
corresponds to varying the baryonic masses of the NSs. Configurations with supercritical baryonic masses result in 
black holes while those with subcritical baryonic masses result in NSs. The system is carried away from the critical 

f — solution by the 1 unstable mode. Even in the absence of a scale in the Einstein field equations, in Type I critical 
phenomena, the time of departure from the critical solution scales with respect to constant powers of the distance of 

£ — the initial data from the critical threshold. 

As critical phenomena is found in the above study when the polytropic index of the EOS is tuned, it is conjectured 
that NS-like objects undergoing a slow cooling process, accretion and angular momentum loss may exhibit critical 
phenomena. In order to clarify the astrophysical implications this may entail, we extract the relevant time scales from 
the critical solution. Since exact critical solutions in principle stay on the critical point as time proceeds to infinity, 
we perform a perturbative mode analysis and compare the mode frequencies of the critical solution with frequencies 
of specifically non-radiating modes for TOV configurations with the same baryonic mass. 

We next determine the universality of these solutions with respect to a 1-parameter family of initial data, by 
constructing packets of matter using the polytropic EOS (p = Kp T ) with K = 80 and T — 2.0, whose distributions 
are characterized by a Gaussian rotated about the axis of axisymmetry. We perform collision simulations similar to 
that done in [3] using these packets of matter and observe that the critical solution found there constitutes a universal 
attractor. We suggest strong evidence that the critical point found in Jin et al is a limit cycle rather than a fixed 
point and present phase diagrams using appropriate evolution variables in the 3+1 split. All parameters and variables 
presented in this work are in G = c = 1 units. 

Sec. 2. Perturbative Mode Analysis. In this analysis, we first obtain the time scale of attraction of the system 
toward the critical solution and the dominant frequencies of the oscillations during the period of the critical solution. 
Fig. 1 shows the evolution of the lapse, density and 4-scalar curvature at the center of collision. We first obtain, 
via a Fourier analysis, the dominant frequencies of the oscillations of these 3 evolution variables. Fig. 2 indicates 2 
dominant frequencies, namely, w% = O.15M0 and 0J2 — 0.23Mq, which correspond to oscillation periods T\ — 0.21tos 
and T 2 = 0.13ms respectively, that are universal across all 3 variables. Using the dominant frequencies, we then fit 
the entire critical solution with the following equations: 
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where a, b, c and d are fit constants, a^, p^, and R x are respectively the averages of the lapse, density and 4-scalar 
curvature throughout the oscillation phase of the critical solution, ratio a = 0.2, ratio p = ration = 0-4 are the ratio of 
the amplitudes of lu 2 over io\ obtained from the Fourier transforms of the respective evolution variables. The fit results 
are shown in Fig. 3. From the fit, we find that A — 0.09, which corresponds to the attraction time, t = 0.05ms, of the 
system toward the critical solution, is universal across the 3 evolution variables. We further observe that the attraction 
time is similar to the departure time found in |2J . This indicates the existence of a complex pair of eigenvalues of the 
growing unstable mode when linearization in a small neighbourhood of the critical solution is performed. We note 
that these time scales are of an order of magnitude smaller than the cooling time scales of NSs commonly reported 
in the astronomy and astrophysics literature. Therefore, there is a high possibility that real astrophysical systems of 
cooling NSs pass through the threshold of critical collapse. Further, systems that undergo slow accretion and angular 
momentum loss with time scales larger than that reported here for the NS critical solution may also pass through the 
threshold and exhibit critical behavior. 

We next perform a perturbative mode analysis on TOV configurations that have the same baryonic mass as the 
object that undergoes critical collapse, namely at baryonic mass Mf, = 1.5393M Q . We note that the exact critical 
solution is non-radiative in principle, and therefore, we focus on the 1=0,1 modes. Using the following perturbation 
metric on a TOV background [I]: 
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with a gauge choice that leaves the exterior solution invariantly spherically symmetric, namely G = K = 0, 
ho = hi = 0, we obtain the following perturbed line elements: 
1=0 mode: 



ds 2 



e u (H - l)dt 2 + e x {H 2 + l)dr 2 + r 2 (d9 2 + : 



: ) 



(5) 




a coefficients 
b coefficients 



a 



R 



0.2 0.4 0.6 0.6 

frequency, CO 



FIG. 2: 



1=1 mode: 
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We further refer to the various past established literature [5], [B], [7] and [5], reconcile the differing conventions and 
confirm their results with our mode analysis. 

For the 1=0 mode, we denote the fluid perturbation as: 



£ = -r- 2 e- x/2 W. 



(7) 



Assuming harmonic dependence on time for the fluid perturbation as the only degree of freedom in this calculation, 
we thus write: 



W(r,t) = W(r)e i 



We then solve the following equations: 
SGI = 87rT° : 



r- l e- x B 2 .t - %-K{p + p)r- z e-\/2W tt = 
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and obtain the following lst-order ODE: 
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where Z = ( £ r)(^ — 2r 1 ), imposing the boundary condition that the pressure perturbation goes to zero at the TOV 
star surface, which, according to [8], is equivalent to: 



(r- 2 e- x ' 2 W)' = ( 7 i?)" 1 (4 + (^r)e x W + uj 2 (^-)e- v )(r- 2 e- x / 2 W), 
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where M is the total ADM mass of the TOV star. From a simple shoot-and-match calculation of a finite-differenced 
TOV configuration as mentioned above, we find that the fundamental frequency is oj\ p = 6.858 x 10 _3 M Q . 
Similarly, for the 1=1 mode, we have: 
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as the fluid perturbations, and write their corresponding harmonic dependences: 

H (r,t) = H a {r)e^ v 
Hi(r,t) = H x (r)e iuJt 
H 2 (r,t) = H 2 (r)e iut 
W(r,t) = W(r)e lu>t _ 

We solve the following equations: 
SG° = 87rT : 

H' 2 + (r _1 (! + eA ) ~ A ' + 4vrre A (p + p))H 2 = 87rre A [(p + p) r - 2 e- x/2 W' + r~ 2 e- x/2 p'W + 2 r - 2 {p + p)V 
5G{=8-kT(: 

S Viw) = 0: 

-{p + p)e' v u 2 V 
to obtain a 3rd-order ODE: 
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We then impose the boundary conditions that the space-time perturbation functions vanish at the TOV star surface 
whilst at the center of the star, the perturbation functions approach a finite value, namely: 
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where w and ft. are arbitrary constants. Via a variational principle calculation following [5], the fundamental frequency 
for this mode is found to be uj 2p — 1.813 x 1O _2 M0. 

We note that ui\ p and to 2p are both of 1 and 2 orders of magnitude smaller than u)\ and lu 2 . Since u>i p and lo 2p 
are the fundamental frequencies of the non-radiating modes of equilibrium configurations, our mode analysis confirms 
that the configuration that undergoes critical collapse is far from equilibrium, and thus a perturbative mode analysis 
has to be performed on non-equilibrium rather than on TOV background space-times in order to determine the nature 
of the oscillation frequencies exhibited by the critical solution. 
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FIG. 4: The Gaussian distributions are constructed so as to maintain the baryonic mass of the system as a constant. Therefore, 
when the Gaussian heights are increased, their widths are correspondingly decreased. Configurations 1 and 2 are less compact 
distributions whereas that of 3 and 4 are more compact. 

Sec. 3. Universality. We note that the initial data constructed by varying the boost velocity, initial separation 
and polytropic EOS are very close to each other and are within the perturbative region of the NS critical solution 
found in the former study. The universality of the critical solution with respect to a 1-parameter family of initial data 
is thus considered unclear [5], Therefore, in this section, we present a completely different family of initial data where 
the matter field consists of packets of matter whose densities are characterized by Gaussian distributions rotated 
about the axis of axisymmetry and a polytropic EOS (p = Kp r ) with K — 80 and V — 2.0, where the tuning of the 
height and width of the Gaussian density distributions produce a 1-parameter family of initial data. Using the same 
numerical setup as in [2], namely T- freezing shift and 1+log slicing 3+1 BSSN evolution on a similar grid, 2 Gaussian 
packets are boosted to a head-on collision and the behavior of the axisymmetric object formed is observed. Although 
the system is in a non-equilibrium state at the initial time, we have checked that their instability time scale, which is 
described by the oscillation time scale of a single Gaussian packet about an equilibrium TOV configuration, is much 
larger than the time scale of their merging and collapse, that the effect of the instability can be considered negligible. 
Table 1 and Fig. 4 show the configurations that are constructed and their density distributions along the direction 



Configuration Mb 



d 



NS 1.5 0.00056595 27.6 0.15 

1 1.5 0.00038698 29.6 0.1 

2 1.5 0.00039192 29.6 0.1 

3 1.6 0.00055501 27.6 0.12 

4 1.6 0.00064912 27.6 0.12 



TABLE I: Mb is the baryonic mass, p c is the NS central density/height of the Gaussian density distribution at t = 0, d is the 
center-to-center separation between the NSs/Gaussian packets, and v z is the boost velocity of the NSs/Gaussian packets along 
the z-direction of the grid. 



of collision respectively. In Fig. 5, we draw the phase trajectories for these configurations using the density at the 
center of the grid, the trace of the extrinsic curvature at the coordinate (2.1,0.06,0.06) and the z-component of the 
3-momentum, Si, of the fluid element on the edge of the density contour bounded from below by 7.5 x 10 -3 and 
compare them with that for the NS configuration in [2] (the next section contains the reasons why these evolution 
variables are chosen for the phase diagrams). We observe that these configurations are all attracted to the same 
attraction basin in the phase diagrams as that of the NS critical solution. When the Gaussian heights exceed the 
critical values for each of the Gaussian packet configurations, NS-like objects are formed, while heights lower than 
the critical values produce black holes. This indicates the universality of the NS critical solution with respect to a 
1-parameter family of initial data, and that the NS attraction basin constitutes a universal attractor. 

Sec. 4. Phase Space Analysis. In this section, we present a probe on the structure of the critical solution in a 
phase space framework. We note that the initial data configurations as presented in the previous section provide us 
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FIG. 5: The trajectories of all the configurations start on the left part and are attracted to the attracting basin on the right 
part. The arrows point to the differing starting points of the trajectories. K denotes the trace of the extrinsic curvature at the 
coordinate (2. f ,0.06,0.06) and S z denotes the z-component of the 3-momentum, Si, of the fluid element on 7.5 x 10~ 3 density 
contour, respectively. The natural logarithm is taken of the density at the center of the grid, namely lnp c , so as to greater 
distinguish the differing starting points of the trajectories in the phase diagrams. 
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with another degree of freedom in probing the structure of the critical collapse threshold in a more global scale. For 
this purpose, we vary the boost velocity of the Gaussian packets together with the Gaussian heights and widths while 
maintaining the baryonic mass of the system at M& = 1.5393. In Fig. 6, we plot boost velocity with respect to Gaussian 
height. We observe that there is a turning point in this phase diagram at approximately (p c ,v) = (0.00064,0.12), 
which corresponds to the NS critical solution. An initial data configuration to the left of the critical surface collapses 
to a black hole and one to the right results in an NS-like object. This critical surface is an indication of the NS 
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FIG. 7: Left legend shows contour densities at which measures are performed. Polar and equatorial proper radii shown on the 
right are measured at the 5.0 x 10~ 4 density contour. The polar direction is taken to be the z-axis on the xz-plane of the grid 
and the equatorial direction is taken to be the x-axis on the same plane. The equatorial direction is thus along the direction 
of collision of the NSs. 



attraction basin as commented on in Sec. 3. Beyond both edges of the critical surface, no critical collapse behavior 
is observed. We therefore note that configurations with Gaussian packets with compactness outside the range of 
0.00033 < p c < 0.00064 are not attracted to the critical surface at all. However, due to the parabolic structure of this 
critical surface, configurations that are less compact than that with p c ~ 0.00064 are more prone to undergo critical 
collapse than those that are more compact. We therefore note that not only the baryonic mass, but also the mass 
distribution, contributes to the dynamical mechanism that drives critical collapse. 

This is further reinforced by the observation of the presence of a stationary envelope surrounding the core of the 
merged object which undergoes critical collapse. For the NS case, we measure the baryonic mass enclosed within 
different density contours throughout the merged object, the proper polar circumferences of these density contours, 
and their equatorial and proper radii, and observe the pattern of their evolution throughout the critical collapse. 
The measure of the baryonic mass enclosed within different density contours throughout the merged object has the 
advantage of being a fully-geometric measure. In Fig. 7, we observe that this measure exhibit oscillations that are 
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correlated between different regions of the merged object, whereas the corresponding proper polar circumferences 
remain basically stationary throughout the critical collapse phase ie. matter fluxes exist through the density contour 
with density threshold p ~ 3.5 x 1CP 3 . The proper radii however exhibit prolonged oscillations only within the region 
bounded by this specific density contour. The oscillations of these radii outside this region exhibit damping into a 
stationary state. These oscillations are correlated between the polar and equatorial directions ie. the polar proper 
radius is in phase with the oscillation of the central density while the equatorial radius is out of phase. 

We also note that the oscillation characteristic is evident across the NS and Gaussian packet critical collapses. 
Further, the oscillation phases of certain evolution variables (as those mentioned in Sec. 3) correlate with each other 
such that they form closed periodic orbits in the phase diagrams, similar to the circle in the phase diagram of a simple 
pendulum. In infinite-dimensional dynamical systems theory |10j , we note that phase space trajectory flows are tangent 
to invariant sets which can be divided into 3 types, ie. the center, stable and unstable subspaces. Stable manifolds 
contain flows that move toward the center manifold whilst unstable manifolds contain flows that move away. These 
manifolds are tangent to their respective subspaces. As a direct analogy, the critical collapse threshold corresponds to 
a 3c" -1 ) center manifold that is bounded by stable and unstable ffi 1 manifolds. When a parameter of the initial data 
is varied, bifurcation occurs and solutions are carried away from the center manifold via the unstable manifold. Since 
stable manifolds dcfinitionally contain flows that move toward the center manifold, bifurcation dcfinitionally could 
not occur across stable manifolds. In both the NS and Gaussian packet cases, bifurcation is seen to occur across an 
oscillatory threshold. This strongly suggests that the center manifold in these cases contain a limit cycle rather than 
a fixed point. 

Sec. 5. Conclusions. In this work, we confirm that NS-like objects undergoing critical collapse are in a non- 
equilibrium state and therefore, their oscillations cannot be described by a perturbative mode analysis on equilibrium 
background space-times. The oscillation frequencies of these objects are found to be of 1 and 2 orders of magnitude 
smaller than that for equilibrium configurations. We note that the time scales of these oscillations, T\ — 0.21ms and 
T-2 = 0.13ms, as well as the time scale of solutions approaching and leaving the critical collapse threshold, t — 0.05ms, 
are also of 1 order of magnitude smaller than the cooling time scales of NSs. This suggests a high possibility that 
NS systems undergoing a slow softening of the EOS, accretion or angular momentum loss undergo critical collapse. 
Critical phenomena in gravitational collapses of real astrophysical systems of NSs will thus have an effect on the 
emission of gravitational radiation that can be detected by LIGO, LISA etc. Work to investigate critical collapses in 
rotating NS-like objects, which greater resemble real NS systems, will be reported elsewhere by the authors KJJ and 
WMS. 

We also confirm, in this work, that the NS critical solution constitutes a universal attraction basin for 1-parameter 
families of NS-like initial data. The unstable mode carrying solutions away from this critical collapse threshold is 
characterized by the complex pair of eigenvalues ±0.09. We present the critical surface resulting from bifurcations of a 
1-parameter family of NS-like initial data and deduce that less compact mass distributions of NS-like objects are more 
prone to undergo critical collapses than more compact ones. With phase diagrams constructed from geometrically- 
invariant quantities and evolution variables of the 3+1 BSSN formulation using T-freezing shift and 1+log slicing, 
we present strong evidence that the NS critical collapse attractor is a limit cycle rather than a fixed point. Since 
an exact critical solution is by definition non-radiative, the presence of a limit cycle on the NS critical collapse 
threshold together with the observation of critical collapse in rotating NS-like systems, poses an interesting question 
of how a gravitational system maintains oscillations and rotations without emitting gravitational radiation. Further 
investigations of this question will be closely related to perturbative mode analyses on non-equilibrium configurations. 
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